import numpy
import numpy as np
import scipy.io.wavfile
import matplotlib.pyplot as plt
from scipy.fftpack import dct

show_or_not = True

# 预处理start
sample_rate, signal = scipy.io.wavfile.read('./yhs.wav')  # 这个音频有多个声道！！！

signal = signal[0:int(3.5 * sample_rate)]   # 我们只取第一个声道的前3.5s

signal = signal[::int(sample_rate/8000)]  # 不规范地转成采样率为8000
sample_rate = 8000

pre_emphasis = 0.97
emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])

frame_size = 0.025
frame_stride = 0.01
frame_length, frame_step = frame_size * sample_rate, frame_stride * sample_rate  # 从秒转换为采样点
signal_length = len(emphasized_signal)
frame_length = int(round(frame_length))
frame_step = int(round(frame_step))
# 除了第一帧额外的帧
num_frames = int(numpy.ceil(float(numpy.max([signal_length - frame_length, 0])) / frame_step))

pad_signal_length = num_frames * frame_step + frame_length
z = numpy.zeros((pad_signal_length - signal_length))
# 把不足一帧的部分补成一帧
pad_signal = numpy.append(emphasized_signal, z)

# 拿到每一个frame的所有下标
indices = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1)) + numpy.tile(numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T
# 根据下标求出每个frame的各个位置对应数据
frames = pad_signal[indices.astype(numpy.int32, copy=False)]

# 每一个窗口对应位置成窗的值
frames *= numpy.hamming(frame_length)
# frames *= 0.54 - 0.46 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1))  # 内部实现

# 预处理end

# FFT start

NFFT = 512
mag_frames = numpy.absolute(numpy.fft.rfft(frames, NFFT))   # fft的幅度(magnitude) rfft只保留了fft图像的一半 把对称的部分删去了

if show_or_not:
    plt.plot(mag_frames[-1])  # 这里就只看最后一帧
    plt.show()

if show_or_not:
    plt.imshow(mag_frames.T)
    plt.show()

# FFT end

# 功率谱 start

pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2))  # 可以不要 如果要功率谱是这么搞
pow_frames = mag_frames  # 这里不要

# 功率谱 end

nfilt = 40
low_freq_mel = 0
high_freq_mel = (2595 * np.log10(1 + (sample_rate / 2) / 700))  # 将Hz转换为Mel
# 我们要做40个滤波器组，为此需要42个点（毕竟第一个需要f0和最后一个需要f41），这意味着在们需要low_freq_mel和high_freq_mel之间线性间隔40个点
mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # 使得Mel scale间距相等
hz_points = (700 * (10 ** (mel_points / 2595) - 1))  # 将Mel转换回-Hz
bins = np.floor((NFFT + 1) * hz_points / sample_rate)  # 和hz的分布成正比 拿到各个三角滤波的

fbank = np.zeros((nfilt, int(np.floor(NFFT / 2 + 1))))
for m in range(1, nfilt + 1):
    f_m_minus = int(bins[m - 1])  # 左
    f_m = int(bins[m])  # 中
    f_m_plus = int(bins[m + 1])  # 右

    for k in range(f_m_minus, f_m):
        fbank[m - 1, k] = (k - bins[m - 1]) / (bins[m] - bins[m - 1])
    for k in range(f_m, f_m_plus):
        fbank[m - 1, k] = (bins[m + 1] - k) / (bins[m + 1] - bins[m])

filter_banks = np.dot(pow_frames, fbank.T)  # 拿到每个时刻的每个滤波器的结果

if show_or_not:
    for i in fbank:
        plt.plot(i)
    plt.show()

if show_or_not:
    plt.imshow(filter_banks.T)
    plt.show()
